import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import netCDF4
import pandas as pd
import scipy.interpolate as interp
%matplotlib inline
# Colormap selection
xr.set_options(cmap_divergent='bwr', cmap_sequential='turbo')
<xarray.core.options.set_options at 0x7fd80c14d5b0>
mfdataDIR1 = 'data/GPM/2009/3B-MO.MS.MRG.3IMERG.*.V06B.HDF5.SUB.nc4'
mfdataDIR2 = 'data/GPM/2019/3B-MO.MS.MRG.3IMERG.*.V06B.HDF5.SUB.nc4'
ds1 = xr.open_mfdataset(mfdataDIR1, parallel=True)
ds2 = xr.open_mfdataset(mfdataDIR2, parallel=True)
ds1
<xarray.Dataset>
Dimensions: (lat: 1800, lon: 3600, time: 12)
Coordinates:
* time (time) datetime64[ns] 2009-01-01 2009-02-01 ... 2009-12-01
* lon (lon) float32 -179.9 -179.9 -179.8 ... 179.8 179.9 179.9
* lat (lat) float32 -89.95 -89.85 -89.75 ... 89.75 89.85 89.95
Data variables:
precipitation (time, lon, lat) float32 dask.array<chunksize=(1, 3600, 1800), meta=np.ndarray>
Attributes:
CDI: Climate Data Interface version 1....
Conventions: CF-1.6
Original_Producer_Metadata_FileHeader: DOI=10.5067/GPM/IMERG/3B-MONTH/06...
Original_Producer_Metadata_FileInfo: DataFormatVersion=6a;\nTKCodeBuil...
Original_Producer_Metadata_GridHeader: BinMethod=ARITHMETIC_MEAN;\nRegis...
InputPointer: 3B-MO.MS.MRG.3IMERG.20090101-S000...
history_L34RS: 'Created by L34RS v1.4.2 @ NASA G...
CDO: Climate Data Operators version 1....array(['2009-01-01T00:00:00.000000000', '2009-02-01T00:00:00.000000000',
'2009-03-01T00:00:00.000000000', '2009-04-01T00:00:00.000000000',
'2009-05-01T00:00:00.000000000', '2009-06-01T00:00:00.000000000',
'2009-07-01T00:00:00.000000000', '2009-08-01T00:00:00.000000000',
'2009-09-01T00:00:00.000000000', '2009-10-01T00:00:00.000000000',
'2009-11-01T00:00:00.000000000', '2009-12-01T00:00:00.000000000'],
dtype='datetime64[ns]')array([-179.95, -179.85, -179.75, ..., 179.75, 179.85, 179.95],
dtype=float32)array([-89.95, -89.85, -89.75, ..., 89.75, 89.85, 89.95], dtype=float32)
|
# make preciptation rate to preciptation
def convert_to_precipitaion(ds):
temp = ds * 24
# temp = temp.to_dataset()
return temp
ds1 = convert_to_precipitaion(ds1)
# Transpose the data to get lat first and lon after -
ds1 = ds1.transpose("time", "lat", "lon")
ds1_ind = ds1.sel(lat=slice(7,36), lon=slice(67,98)).dropna("time")
# Wrap it into a simple function
def season_mean(ds, calendar='standard'):
# Make a DataArray with the number of days in each month, size = len(time)
month_length = ds.time.dt.days_in_month
# Calculate the weights by grouping by 'time.season'
weights = month_length.groupby('time.season') / month_length.groupby('time.season').sum()
# Test that the sum of the weights for each season is 1.0
np.testing.assert_allclose(weights.groupby('time.season').sum().values, np.ones(4))
# Calculate the weighted average
return (ds * weights).groupby('time.season').sum(dim='time')
# Get seasonal mean
ds1_ind_sm = season_mean(ds1_ind)
ds1_ind_sm
<xarray.Dataset>
Dimensions: (lat: 290, lon: 310, season: 4)
Coordinates:
* lon (lon) float32 67.05 67.15 67.25 67.35 ... 97.75 97.85 97.95
* lat (lat) float32 7.05 7.15 7.25 7.35 ... 35.65 35.75 35.85 35.95
* season (season) object 'DJF' 'JJA' 'MAM' 'SON'
Data variables:
precipitation (season, lat, lon) float64 dask.array<chunksize=(1, 290, 310), meta=np.ndarray>array([67.05, 67.15, 67.25, ..., 97.75, 97.85, 97.95], dtype=float32)
array([ 7.05, 7.15, 7.25, ..., 35.75, 35.85, 35.95], dtype=float32)
array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
|
# Convert to dataarray
da1 = ds1_ind_sm.precipitation
da1
<xarray.DataArray 'precipitation' (season: 4, lat: 290, lon: 310)> dask.array<concatenate, shape=(4, 290, 310), dtype=float64, chunksize=(1, 290, 310), chunktype=numpy.ndarray> Coordinates: * lon (lon) float32 67.05 67.15 67.25 67.35 ... 97.65 97.75 97.85 97.95 * lat (lat) float32 7.05 7.15 7.25 7.35 7.45 ... 35.65 35.75 35.85 35.95 * season (season) object 'DJF' 'JJA' 'MAM' 'SON'
|
array([67.05, 67.15, 67.25, ..., 97.75, 97.85, 97.95], dtype=float32)
array([ 7.05, 7.15, 7.25, ..., 35.75, 35.85, 35.95], dtype=float32)
array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
da1.sel(season = 'JJA').plot()
<matplotlib.collections.QuadMesh at 0x7fd7b44ec370>
import geopandas as gpd
from rasterio import features
from affine import Affine
def transform_from_latlon(lat, lon):
""" input 1D array of lat / lon and output an Affine transformation
"""
lat = np.asarray(lat)
lon = np.asarray(lon)
trans = Affine.translation(lon[0], lat[0])
scale = Affine.scale(lon[1] - lon[0], lat[1] - lat[0])
return trans * scale
def rasterize(shapes, coords, latitude='lat', longitude='lon',
fill=np.nan, **kwargs):
"""Rasterize a list of (geometry, fill_value) tuples onto the given
xray coordinates. This only works for 1d latitude and longitude
arrays.
usage:
-----
1. read shapefile to geopandas.GeoDataFrame
`states = gpd.read_file(shp_dir+shp_file)`
2. encode the different shapefiles that capture those lat-lons as different
numbers i.e. 0.0, 1.0 ... and otherwise np.nan
`shapes = (zip(states.geometry, range(len(states))))`
3. Assign this to a new coord in your original xarray.DataArray
`ds['states'] = rasterize(shapes, ds.coords, longitude='X', latitude='Y')`
arguments:
---------
: **kwargs (dict): passed to `rasterio.rasterize` function
attrs:
-----
:transform (affine.Affine): how to translate from latlon to ...?
:raster (numpy.ndarray): use rasterio.features.rasterize fill the values
outside the .shp file with np.nan
:spatial_coords (dict): dictionary of {"X":xr.DataArray, "Y":xr.DataArray()}
with "X", "Y" as keys, and xr.DataArray as values
returns:
-------
:(xr.DataArray): DataArray with `values` of nan for points outside shapefile
and coords `Y` = latitude, 'X' = longitude.
"""
transform = transform_from_latlon(coords['lat'], coords['lon'])
out_shape = (len(coords['lat']), len(coords['lon']))
raster = features.rasterize(shapes, out_shape=out_shape,
fill=fill, transform=transform,
dtype=float, **kwargs)
spatial_coords = {latitude: coords['lat'], longitude: coords['lon']}
return xr.DataArray(raster, coords=spatial_coords, dims=('lat', 'lon'))
def add_shape_coord_from_data_array(xr_da, shp_path, coord_name):
""" Create a new coord for the xr_da indicating whether or not it
is inside the shapefile
Creates a new coord - "coord_name" which will have integer values
used to subset xr_da for plotting / analysis/
Usage:
-----
precip_da = add_shape_coord_from_data_array(precip_da, "awash.shp", "awash")
awash_da = precip_da.where(precip_da.awash==0, other=np.nan)
"""
# 1. read in shapefile
shp_gpd = gpd.read_file(shp_path)
# 2. create a list of tuples (shapely.geometry, id)
# this allows for many different polygons within a .shp file (e.g. States of US)
shapes = [(shape, n) for n, shape in enumerate(shp_gpd.geometry)]
# 3. create a new coord in the xr_da which will be set to the id in `shapes`
xr_da[coord_name] = rasterize(shapes, xr_da.coords,
longitude='longitude', latitude='latitude')
return xr_da
shp_dir = './shapefiles/'
da1_ind = add_shape_coord_from_data_array(da1, shp_dir, "awash")
awash_da1 = da1_ind.where(da1_ind.awash==0, other=np.nan)
awash_da1.sel(season="JJA").plot()
<matplotlib.collections.QuadMesh at 0x7fd782e37220>
awash_da1
<xarray.DataArray 'precipitation' (season: 4, lat: 290, lon: 310)>
dask.array<where, shape=(4, 290, 310), dtype=float64, chunksize=(1, 290, 310), chunktype=numpy.ndarray>
Coordinates:
* lon (lon) float32 67.05 67.15 67.25 67.35 ... 97.65 97.75 97.85 97.95
* lat (lat) float32 7.05 7.15 7.25 7.35 ... 35.65 35.75 35.85 35.95
* season (season) object 'DJF' 'JJA' 'MAM' 'SON'
latitude (lat) float32 7.05 7.15 7.25 7.35 ... 35.65 35.75 35.85 35.95
longitude (lon) float32 67.05 67.15 67.25 67.35 ... 97.65 97.75 97.85 97.95
awash (lat, lon) float64 nan nan nan nan nan ... nan nan nan nan nan
|
array([67.05, 67.15, 67.25, ..., 97.75, 97.85, 97.95], dtype=float32)
array([ 7.05, 7.15, 7.25, ..., 35.75, 35.85, 35.95], dtype=float32)
array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
array([ 7.05, 7.15, 7.25, 7.35, 7.45, 7.55, 7.65, 7.75, 7.85,
7.95, 8.05, 8.15, 8.25, 8.35, 8.45, 8.55, 8.65, 8.75,
8.85, 8.95, 9.05, 9.15, 9.25, 9.35, 9.45, 9.55, 9.65,
9.75, 9.85, 9.95, 10.05, 10.15, 10.25, 10.35, 10.45, 10.55,
10.65, 10.75, 10.85, 10.95, 11.05, 11.15, 11.25, 11.35, 11.45,
11.55, 11.65, 11.75, 11.85, 11.95, 12.05, 12.15, 12.25, 12.35,
12.45, 12.55, 12.65, 12.75, 12.85, 12.95, 13.05, 13.15, 13.25,
13.35, 13.45, 13.55, 13.65, 13.75, 13.85, 13.95, 14.05, 14.15,
14.25, 14.35, 14.45, 14.55, 14.65, 14.75, 14.85, 14.95, 15.05,
15.15, 15.25, 15.35, 15.45, 15.55, 15.65, 15.75, 15.85, 15.95,
16.05, 16.15, 16.25, 16.35, 16.45, 16.55, 16.65, 16.75, 16.85,
16.95, 17.05, 17.15, 17.25, 17.35, 17.45, 17.55, 17.65, 17.75,
17.85, 17.95, 18.05, 18.15, 18.25, 18.35, 18.45, 18.55, 18.65,
18.75, 18.85, 18.95, 19.05, 19.15, 19.25, 19.35, 19.45, 19.55,
19.65, 19.75, 19.85, 19.95, 20.05, 20.15, 20.25, 20.35, 20.45,
20.55, 20.65, 20.75, 20.85, 20.95, 21.05, 21.15, 21.25, 21.35,
21.45, 21.55, 21.65, 21.75, 21.85, 21.95, 22.05, 22.15, 22.25,
22.35, 22.45, 22.55, 22.65, 22.75, 22.85, 22.95, 23.05, 23.15,
23.25, 23.35, 23.45, 23.55, 23.65, 23.75, 23.85, 23.95, 24.05,
24.15, 24.25, 24.35, 24.45, 24.55, 24.65, 24.75, 24.85, 24.95,
25.05, 25.15, 25.25, 25.35, 25.45, 25.55, 25.65, 25.75, 25.85,
25.95, 26.05, 26.15, 26.25, 26.35, 26.45, 26.55, 26.65, 26.75,
26.85, 26.95, 27.05, 27.15, 27.25, 27.35, 27.45, 27.55, 27.65,
27.75, 27.85, 27.95, 28.05, 28.15, 28.25, 28.35, 28.45, 28.55,
28.65, 28.75, 28.85, 28.95, 29.05, 29.15, 29.25, 29.35, 29.45,
29.55, 29.65, 29.75, 29.85, 29.95, 30.05, 30.15, 30.25, 30.35,
30.45, 30.55, 30.65, 30.75, 30.85, 30.95, 31.05, 31.15, 31.25,
31.35, 31.45, 31.55, 31.65, 31.75, 31.85, 31.95, 32.05, 32.15,
32.25, 32.35, 32.45, 32.55, 32.65, 32.75, 32.85, 32.95, 33.05,
33.15, 33.25, 33.35, 33.45, 33.55, 33.65, 33.75, 33.85, 33.95,
34.05, 34.15, 34.25, 34.35, 34.45, 34.55, 34.65, 34.75, 34.85,
34.95, 35.05, 35.15, 35.25, 35.35, 35.45, 35.55, 35.65, 35.75,
35.85, 35.95], dtype=float32)array([67.05, 67.15, 67.25, 67.35, 67.45, 67.55, 67.65, 67.75, 67.85,
67.95, 68.05, 68.15, 68.25, 68.35, 68.45, 68.55, 68.65, 68.75,
68.85, 68.95, 69.05, 69.15, 69.25, 69.35, 69.45, 69.55, 69.65,
69.75, 69.85, 69.95, 70.05, 70.15, 70.25, 70.35, 70.45, 70.55,
70.65, 70.75, 70.85, 70.95, 71.05, 71.15, 71.25, 71.35, 71.45,
71.55, 71.65, 71.75, 71.85, 71.95, 72.05, 72.15, 72.25, 72.35,
72.45, 72.55, 72.65, 72.75, 72.85, 72.95, 73.05, 73.15, 73.25,
73.35, 73.45, 73.55, 73.65, 73.75, 73.85, 73.95, 74.05, 74.15,
74.25, 74.35, 74.45, 74.55, 74.65, 74.75, 74.85, 74.95, 75.05,
75.15, 75.25, 75.35, 75.45, 75.55, 75.65, 75.75, 75.85, 75.95,
76.05, 76.15, 76.25, 76.35, 76.45, 76.55, 76.65, 76.75, 76.85,
76.95, 77.05, 77.15, 77.25, 77.35, 77.45, 77.55, 77.65, 77.75,
77.85, 77.95, 78.05, 78.15, 78.25, 78.35, 78.45, 78.55, 78.65,
78.75, 78.85, 78.95, 79.05, 79.15, 79.25, 79.35, 79.45, 79.55,
79.65, 79.75, 79.85, 79.95, 80.05, 80.15, 80.25, 80.35, 80.45,
80.55, 80.65, 80.75, 80.85, 80.95, 81.05, 81.15, 81.25, 81.35,
81.45, 81.55, 81.65, 81.75, 81.85, 81.95, 82.05, 82.15, 82.25,
82.35, 82.45, 82.55, 82.65, 82.75, 82.85, 82.95, 83.05, 83.15,
83.25, 83.35, 83.45, 83.55, 83.65, 83.75, 83.85, 83.95, 84.05,
84.15, 84.25, 84.35, 84.45, 84.55, 84.65, 84.75, 84.85, 84.95,
85.05, 85.15, 85.25, 85.35, 85.45, 85.55, 85.65, 85.75, 85.85,
85.95, 86.05, 86.15, 86.25, 86.35, 86.45, 86.55, 86.65, 86.75,
86.85, 86.95, 87.05, 87.15, 87.25, 87.35, 87.45, 87.55, 87.65,
87.75, 87.85, 87.95, 88.05, 88.15, 88.25, 88.35, 88.45, 88.55,
88.65, 88.75, 88.85, 88.95, 89.05, 89.15, 89.25, 89.35, 89.45,
89.55, 89.65, 89.75, 89.85, 89.95, 90.05, 90.15, 90.25, 90.35,
90.45, 90.55, 90.65, 90.75, 90.85, 90.95, 91.05, 91.15, 91.25,
91.35, 91.45, 91.55, 91.65, 91.75, 91.85, 91.95, 92.05, 92.15,
92.25, 92.35, 92.45, 92.55, 92.65, 92.75, 92.85, 92.95, 93.05,
93.15, 93.25, 93.35, 93.45, 93.55, 93.65, 93.75, 93.85, 93.95,
94.05, 94.15, 94.25, 94.35, 94.45, 94.55, 94.65, 94.75, 94.85,
94.95, 95.05, 95.15, 95.25, 95.35, 95.45, 95.55, 95.65, 95.75,
95.85, 95.95, 96.05, 96.15, 96.25, 96.35, 96.45, 96.55, 96.65,
96.75, 96.85, 96.95, 97.05, 97.15, 97.25, 97.35, 97.45, 97.55,
97.65, 97.75, 97.85, 97.95], dtype=float32)array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]])# Plotting
fig = plt.figure(figsize=(20, 15))
fig.tight_layout()
titles = ["MAM", "JJA", "SON", "DJF"]
for i,season in enumerate(titles):
ax = fig.add_subplot(2, 2, i+1, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
awash_da1.sel(season=titles[i]).plot()
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("Months"+ " " + "-" + " " + "("+titles[i]+")", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
fig.suptitle('Precipitation over India (in mm) year 2009', fontsize=20, y=0.95)
plt.savefig('./images/GPM2009.png')
ds2
<xarray.Dataset>
Dimensions: (lat: 1800, lon: 3600, time: 12)
Coordinates:
* time (time) datetime64[ns] 2019-01-01 2019-02-01 ... 2019-12-01
* lon (lon) float32 -179.9 -179.9 -179.8 ... 179.8 179.9 179.9
* lat (lat) float32 -89.95 -89.85 -89.75 ... 89.75 89.85 89.95
Data variables:
precipitation (time, lon, lat) float32 dask.array<chunksize=(1, 3600, 1800), meta=np.ndarray>
Attributes:
CDI: Climate Data Interface version 1....
Conventions: CF-1.6
Original_Producer_Metadata_FileHeader: DOI=10.5067/GPM/IMERG/3B-MONTH/06...
Original_Producer_Metadata_FileInfo: DataFormatVersion=6a;\nTKCodeBuil...
Original_Producer_Metadata_GridHeader: BinMethod=ARITHMETIC_MEAN;\nRegis...
InputPointer: 3B-MO.MS.MRG.3IMERG.20190101-S000...
history_L34RS: 'Created by L34RS v1.4.2 @ NASA G...
CDO: Climate Data Operators version 1....array(['2019-01-01T00:00:00.000000000', '2019-02-01T00:00:00.000000000',
'2019-03-01T00:00:00.000000000', '2019-04-01T00:00:00.000000000',
'2019-05-01T00:00:00.000000000', '2019-06-01T00:00:00.000000000',
'2019-07-01T00:00:00.000000000', '2019-08-01T00:00:00.000000000',
'2019-09-01T00:00:00.000000000', '2019-10-01T00:00:00.000000000',
'2019-11-01T00:00:00.000000000', '2019-12-01T00:00:00.000000000'],
dtype='datetime64[ns]')array([-179.95, -179.85, -179.75, ..., 179.75, 179.85, 179.95],
dtype=float32)array([-89.95, -89.85, -89.75, ..., 89.75, 89.85, 89.95], dtype=float32)
|
# make preciptation rate to preciptation
ds2 = convert_to_precipitaion(ds2)
# Transpose the data to get lat first and lon after -
ds2 = ds2.transpose("time", "lat", "lon")
ds2_ind = ds2.sel(lat=slice(7,36), lon=slice(67,98)).dropna("time")
# Get seasonal mean
ds2_ind_sm = season_mean(ds2_ind)
ds2_ind_sm
<xarray.Dataset>
Dimensions: (lat: 290, lon: 310, season: 4)
Coordinates:
* lon (lon) float32 67.05 67.15 67.25 67.35 ... 97.75 97.85 97.95
* lat (lat) float32 7.05 7.15 7.25 7.35 ... 35.65 35.75 35.85 35.95
* season (season) object 'DJF' 'JJA' 'MAM' 'SON'
Data variables:
precipitation (season, lat, lon) float64 dask.array<chunksize=(1, 290, 310), meta=np.ndarray>array([67.05, 67.15, 67.25, ..., 97.75, 97.85, 97.95], dtype=float32)
array([ 7.05, 7.15, 7.25, ..., 35.75, 35.85, 35.95], dtype=float32)
array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
|
# Convert to dataarray
da2 = ds2_ind_sm.precipitation
da2
<xarray.DataArray 'precipitation' (season: 4, lat: 290, lon: 310)> dask.array<concatenate, shape=(4, 290, 310), dtype=float64, chunksize=(1, 290, 310), chunktype=numpy.ndarray> Coordinates: * lon (lon) float32 67.05 67.15 67.25 67.35 ... 97.65 97.75 97.85 97.95 * lat (lat) float32 7.05 7.15 7.25 7.35 7.45 ... 35.65 35.75 35.85 35.95 * season (season) object 'DJF' 'JJA' 'MAM' 'SON'
|
array([67.05, 67.15, 67.25, ..., 97.75, 97.85, 97.95], dtype=float32)
array([ 7.05, 7.15, 7.25, ..., 35.75, 35.85, 35.95], dtype=float32)
array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
da2.sel(season = 'JJA').plot()
<matplotlib.collections.QuadMesh at 0x7fd783da56a0>
# Masking the data
da2_ind = add_shape_coord_from_data_array(da2, shp_dir, "awash")
awash_da2 = da2_ind.where(da2_ind.awash==0, other=np.nan)
awash_da2.sel(season="JJA").plot()
<matplotlib.collections.QuadMesh at 0x7fd78307e520>
awash_da2
<xarray.DataArray 'precipitation' (season: 4, lat: 290, lon: 310)>
dask.array<where, shape=(4, 290, 310), dtype=float64, chunksize=(1, 290, 310), chunktype=numpy.ndarray>
Coordinates:
* lon (lon) float32 67.05 67.15 67.25 67.35 ... 97.65 97.75 97.85 97.95
* lat (lat) float32 7.05 7.15 7.25 7.35 ... 35.65 35.75 35.85 35.95
* season (season) object 'DJF' 'JJA' 'MAM' 'SON'
latitude (lat) float32 7.05 7.15 7.25 7.35 ... 35.65 35.75 35.85 35.95
longitude (lon) float32 67.05 67.15 67.25 67.35 ... 97.65 97.75 97.85 97.95
awash (lat, lon) float64 nan nan nan nan nan ... nan nan nan nan nan
|
array([67.05, 67.15, 67.25, ..., 97.75, 97.85, 97.95], dtype=float32)
array([ 7.05, 7.15, 7.25, ..., 35.75, 35.85, 35.95], dtype=float32)
array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
array([ 7.05, 7.15, 7.25, 7.35, 7.45, 7.55, 7.65, 7.75, 7.85,
7.95, 8.05, 8.15, 8.25, 8.35, 8.45, 8.55, 8.65, 8.75,
8.85, 8.95, 9.05, 9.15, 9.25, 9.35, 9.45, 9.55, 9.65,
9.75, 9.85, 9.95, 10.05, 10.15, 10.25, 10.35, 10.45, 10.55,
10.65, 10.75, 10.85, 10.95, 11.05, 11.15, 11.25, 11.35, 11.45,
11.55, 11.65, 11.75, 11.85, 11.95, 12.05, 12.15, 12.25, 12.35,
12.45, 12.55, 12.65, 12.75, 12.85, 12.95, 13.05, 13.15, 13.25,
13.35, 13.45, 13.55, 13.65, 13.75, 13.85, 13.95, 14.05, 14.15,
14.25, 14.35, 14.45, 14.55, 14.65, 14.75, 14.85, 14.95, 15.05,
15.15, 15.25, 15.35, 15.45, 15.55, 15.65, 15.75, 15.85, 15.95,
16.05, 16.15, 16.25, 16.35, 16.45, 16.55, 16.65, 16.75, 16.85,
16.95, 17.05, 17.15, 17.25, 17.35, 17.45, 17.55, 17.65, 17.75,
17.85, 17.95, 18.05, 18.15, 18.25, 18.35, 18.45, 18.55, 18.65,
18.75, 18.85, 18.95, 19.05, 19.15, 19.25, 19.35, 19.45, 19.55,
19.65, 19.75, 19.85, 19.95, 20.05, 20.15, 20.25, 20.35, 20.45,
20.55, 20.65, 20.75, 20.85, 20.95, 21.05, 21.15, 21.25, 21.35,
21.45, 21.55, 21.65, 21.75, 21.85, 21.95, 22.05, 22.15, 22.25,
22.35, 22.45, 22.55, 22.65, 22.75, 22.85, 22.95, 23.05, 23.15,
23.25, 23.35, 23.45, 23.55, 23.65, 23.75, 23.85, 23.95, 24.05,
24.15, 24.25, 24.35, 24.45, 24.55, 24.65, 24.75, 24.85, 24.95,
25.05, 25.15, 25.25, 25.35, 25.45, 25.55, 25.65, 25.75, 25.85,
25.95, 26.05, 26.15, 26.25, 26.35, 26.45, 26.55, 26.65, 26.75,
26.85, 26.95, 27.05, 27.15, 27.25, 27.35, 27.45, 27.55, 27.65,
27.75, 27.85, 27.95, 28.05, 28.15, 28.25, 28.35, 28.45, 28.55,
28.65, 28.75, 28.85, 28.95, 29.05, 29.15, 29.25, 29.35, 29.45,
29.55, 29.65, 29.75, 29.85, 29.95, 30.05, 30.15, 30.25, 30.35,
30.45, 30.55, 30.65, 30.75, 30.85, 30.95, 31.05, 31.15, 31.25,
31.35, 31.45, 31.55, 31.65, 31.75, 31.85, 31.95, 32.05, 32.15,
32.25, 32.35, 32.45, 32.55, 32.65, 32.75, 32.85, 32.95, 33.05,
33.15, 33.25, 33.35, 33.45, 33.55, 33.65, 33.75, 33.85, 33.95,
34.05, 34.15, 34.25, 34.35, 34.45, 34.55, 34.65, 34.75, 34.85,
34.95, 35.05, 35.15, 35.25, 35.35, 35.45, 35.55, 35.65, 35.75,
35.85, 35.95], dtype=float32)array([67.05, 67.15, 67.25, 67.35, 67.45, 67.55, 67.65, 67.75, 67.85,
67.95, 68.05, 68.15, 68.25, 68.35, 68.45, 68.55, 68.65, 68.75,
68.85, 68.95, 69.05, 69.15, 69.25, 69.35, 69.45, 69.55, 69.65,
69.75, 69.85, 69.95, 70.05, 70.15, 70.25, 70.35, 70.45, 70.55,
70.65, 70.75, 70.85, 70.95, 71.05, 71.15, 71.25, 71.35, 71.45,
71.55, 71.65, 71.75, 71.85, 71.95, 72.05, 72.15, 72.25, 72.35,
72.45, 72.55, 72.65, 72.75, 72.85, 72.95, 73.05, 73.15, 73.25,
73.35, 73.45, 73.55, 73.65, 73.75, 73.85, 73.95, 74.05, 74.15,
74.25, 74.35, 74.45, 74.55, 74.65, 74.75, 74.85, 74.95, 75.05,
75.15, 75.25, 75.35, 75.45, 75.55, 75.65, 75.75, 75.85, 75.95,
76.05, 76.15, 76.25, 76.35, 76.45, 76.55, 76.65, 76.75, 76.85,
76.95, 77.05, 77.15, 77.25, 77.35, 77.45, 77.55, 77.65, 77.75,
77.85, 77.95, 78.05, 78.15, 78.25, 78.35, 78.45, 78.55, 78.65,
78.75, 78.85, 78.95, 79.05, 79.15, 79.25, 79.35, 79.45, 79.55,
79.65, 79.75, 79.85, 79.95, 80.05, 80.15, 80.25, 80.35, 80.45,
80.55, 80.65, 80.75, 80.85, 80.95, 81.05, 81.15, 81.25, 81.35,
81.45, 81.55, 81.65, 81.75, 81.85, 81.95, 82.05, 82.15, 82.25,
82.35, 82.45, 82.55, 82.65, 82.75, 82.85, 82.95, 83.05, 83.15,
83.25, 83.35, 83.45, 83.55, 83.65, 83.75, 83.85, 83.95, 84.05,
84.15, 84.25, 84.35, 84.45, 84.55, 84.65, 84.75, 84.85, 84.95,
85.05, 85.15, 85.25, 85.35, 85.45, 85.55, 85.65, 85.75, 85.85,
85.95, 86.05, 86.15, 86.25, 86.35, 86.45, 86.55, 86.65, 86.75,
86.85, 86.95, 87.05, 87.15, 87.25, 87.35, 87.45, 87.55, 87.65,
87.75, 87.85, 87.95, 88.05, 88.15, 88.25, 88.35, 88.45, 88.55,
88.65, 88.75, 88.85, 88.95, 89.05, 89.15, 89.25, 89.35, 89.45,
89.55, 89.65, 89.75, 89.85, 89.95, 90.05, 90.15, 90.25, 90.35,
90.45, 90.55, 90.65, 90.75, 90.85, 90.95, 91.05, 91.15, 91.25,
91.35, 91.45, 91.55, 91.65, 91.75, 91.85, 91.95, 92.05, 92.15,
92.25, 92.35, 92.45, 92.55, 92.65, 92.75, 92.85, 92.95, 93.05,
93.15, 93.25, 93.35, 93.45, 93.55, 93.65, 93.75, 93.85, 93.95,
94.05, 94.15, 94.25, 94.35, 94.45, 94.55, 94.65, 94.75, 94.85,
94.95, 95.05, 95.15, 95.25, 95.35, 95.45, 95.55, 95.65, 95.75,
95.85, 95.95, 96.05, 96.15, 96.25, 96.35, 96.45, 96.55, 96.65,
96.75, 96.85, 96.95, 97.05, 97.15, 97.25, 97.35, 97.45, 97.55,
97.65, 97.75, 97.85, 97.95], dtype=float32)array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]])# Plotting
fig = plt.figure(figsize=(20, 15))
fig.tight_layout()
titles = ["MAM", "JJA", "SON", "DJF"]
for i,season in enumerate(titles):
ax = fig.add_subplot(2, 2, i+1, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
awash_da2.sel(season=titles[i]).plot()
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("Months"+ " " + "-" + " " + "("+titles[i]+")", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
fig.suptitle('Precipitation over India (in mm) year 2019', fontsize=20, y=0.95)
plt.savefig('./images/GPM2019.png')
data3 = 'data/IMD/_Clim_Pred_LRF_New_RF25_IMD0p252009.nc'
data4 = 'data/IMD/_Clim_Pred_LRF_New_RF25_IMD0p252019.nc'
ds3 = xr.open_dataset(data3)
ds4 = xr.open_dataset(data4)
ds3
<xarray.Dataset>
Dimensions: (LATITUDE: 129, LONGITUDE: 135, TIME: 365)
Coordinates:
* LONGITUDE (LONGITUDE) float64 66.5 66.75 67.0 67.25 ... 99.5 99.75 100.0
* LATITUDE (LATITUDE) float64 6.5 6.75 7.0 7.25 ... 37.75 38.0 38.25 38.5
* TIME (TIME) datetime64[ns] 2009-01-01 2009-01-02 ... 2009-12-31
Data variables:
RAINFALL (TIME, LATITUDE, LONGITUDE) float32 ...
Attributes:
history: FERRET V6.9 13-Jan-21
Conventions: CF-1.0array([ 66.5 , 66.75, 67. , 67.25, 67.5 , 67.75, 68. , 68.25, 68.5 ,
68.75, 69. , 69.25, 69.5 , 69.75, 70. , 70.25, 70.5 , 70.75,
71. , 71.25, 71.5 , 71.75, 72. , 72.25, 72.5 , 72.75, 73. ,
73.25, 73.5 , 73.75, 74. , 74.25, 74.5 , 74.75, 75. , 75.25,
75.5 , 75.75, 76. , 76.25, 76.5 , 76.75, 77. , 77.25, 77.5 ,
77.75, 78. , 78.25, 78.5 , 78.75, 79. , 79.25, 79.5 , 79.75,
80. , 80.25, 80.5 , 80.75, 81. , 81.25, 81.5 , 81.75, 82. ,
82.25, 82.5 , 82.75, 83. , 83.25, 83.5 , 83.75, 84. , 84.25,
84.5 , 84.75, 85. , 85.25, 85.5 , 85.75, 86. , 86.25, 86.5 ,
86.75, 87. , 87.25, 87.5 , 87.75, 88. , 88.25, 88.5 , 88.75,
89. , 89.25, 89.5 , 89.75, 90. , 90.25, 90.5 , 90.75, 91. ,
91.25, 91.5 , 91.75, 92. , 92.25, 92.5 , 92.75, 93. , 93.25,
93.5 , 93.75, 94. , 94.25, 94.5 , 94.75, 95. , 95.25, 95.5 ,
95.75, 96. , 96.25, 96.5 , 96.75, 97. , 97.25, 97.5 , 97.75,
98. , 98.25, 98.5 , 98.75, 99. , 99.25, 99.5 , 99.75, 100. ])array([ 6.5 , 6.75, 7. , 7.25, 7.5 , 7.75, 8. , 8.25, 8.5 , 8.75,
9. , 9.25, 9.5 , 9.75, 10. , 10.25, 10.5 , 10.75, 11. , 11.25,
11.5 , 11.75, 12. , 12.25, 12.5 , 12.75, 13. , 13.25, 13.5 , 13.75,
14. , 14.25, 14.5 , 14.75, 15. , 15.25, 15.5 , 15.75, 16. , 16.25,
16.5 , 16.75, 17. , 17.25, 17.5 , 17.75, 18. , 18.25, 18.5 , 18.75,
19. , 19.25, 19.5 , 19.75, 20. , 20.25, 20.5 , 20.75, 21. , 21.25,
21.5 , 21.75, 22. , 22.25, 22.5 , 22.75, 23. , 23.25, 23.5 , 23.75,
24. , 24.25, 24.5 , 24.75, 25. , 25.25, 25.5 , 25.75, 26. , 26.25,
26.5 , 26.75, 27. , 27.25, 27.5 , 27.75, 28. , 28.25, 28.5 , 28.75,
29. , 29.25, 29.5 , 29.75, 30. , 30.25, 30.5 , 30.75, 31. , 31.25,
31.5 , 31.75, 32. , 32.25, 32.5 , 32.75, 33. , 33.25, 33.5 , 33.75,
34. , 34.25, 34.5 , 34.75, 35. , 35.25, 35.5 , 35.75, 36. , 36.25,
36.5 , 36.75, 37. , 37.25, 37.5 , 37.75, 38. , 38.25, 38.5 ])array(['2009-01-01T00:00:00.000000000', '2009-01-02T00:00:00.000000000',
'2009-01-03T00:00:00.000000000', ..., '2009-12-29T00:00:00.000000000',
'2009-12-30T00:00:00.000000000', '2009-12-31T00:00:00.000000000'],
dtype='datetime64[ns]')[6356475 values with dtype=float32]
# rename dimension names
ds3_ind = ds3.rename({"LONGITUDE":"lon", "LATITUDE":"lat","TIME":"time"})
ds4_ind = ds4.rename({"LONGITUDE":"lon", "LATITUDE":"lat","TIME":"time"})
# Getting seasonal mean for IMD data
ds3_ind_sm = season_mean(ds3_ind)
ds4_ind_sm = season_mean(ds4_ind)
ds3_ind_sm
<xarray.Dataset>
Dimensions: (lat: 129, lon: 135, season: 4)
Coordinates:
* lon (lon) float64 66.5 66.75 67.0 67.25 ... 99.25 99.5 99.75 100.0
* lat (lat) float64 6.5 6.75 7.0 7.25 7.5 ... 37.5 37.75 38.0 38.25 38.5
* season (season) object 'DJF' 'JJA' 'MAM' 'SON'
Data variables:
RAINFALL (season, lat, lon) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0array([ 66.5 , 66.75, 67. , 67.25, 67.5 , 67.75, 68. , 68.25, 68.5 ,
68.75, 69. , 69.25, 69.5 , 69.75, 70. , 70.25, 70.5 , 70.75,
71. , 71.25, 71.5 , 71.75, 72. , 72.25, 72.5 , 72.75, 73. ,
73.25, 73.5 , 73.75, 74. , 74.25, 74.5 , 74.75, 75. , 75.25,
75.5 , 75.75, 76. , 76.25, 76.5 , 76.75, 77. , 77.25, 77.5 ,
77.75, 78. , 78.25, 78.5 , 78.75, 79. , 79.25, 79.5 , 79.75,
80. , 80.25, 80.5 , 80.75, 81. , 81.25, 81.5 , 81.75, 82. ,
82.25, 82.5 , 82.75, 83. , 83.25, 83.5 , 83.75, 84. , 84.25,
84.5 , 84.75, 85. , 85.25, 85.5 , 85.75, 86. , 86.25, 86.5 ,
86.75, 87. , 87.25, 87.5 , 87.75, 88. , 88.25, 88.5 , 88.75,
89. , 89.25, 89.5 , 89.75, 90. , 90.25, 90.5 , 90.75, 91. ,
91.25, 91.5 , 91.75, 92. , 92.25, 92.5 , 92.75, 93. , 93.25,
93.5 , 93.75, 94. , 94.25, 94.5 , 94.75, 95. , 95.25, 95.5 ,
95.75, 96. , 96.25, 96.5 , 96.75, 97. , 97.25, 97.5 , 97.75,
98. , 98.25, 98.5 , 98.75, 99. , 99.25, 99.5 , 99.75, 100. ])array([ 6.5 , 6.75, 7. , 7.25, 7.5 , 7.75, 8. , 8.25, 8.5 , 8.75,
9. , 9.25, 9.5 , 9.75, 10. , 10.25, 10.5 , 10.75, 11. , 11.25,
11.5 , 11.75, 12. , 12.25, 12.5 , 12.75, 13. , 13.25, 13.5 , 13.75,
14. , 14.25, 14.5 , 14.75, 15. , 15.25, 15.5 , 15.75, 16. , 16.25,
16.5 , 16.75, 17. , 17.25, 17.5 , 17.75, 18. , 18.25, 18.5 , 18.75,
19. , 19.25, 19.5 , 19.75, 20. , 20.25, 20.5 , 20.75, 21. , 21.25,
21.5 , 21.75, 22. , 22.25, 22.5 , 22.75, 23. , 23.25, 23.5 , 23.75,
24. , 24.25, 24.5 , 24.75, 25. , 25.25, 25.5 , 25.75, 26. , 26.25,
26.5 , 26.75, 27. , 27.25, 27.5 , 27.75, 28. , 28.25, 28.5 , 28.75,
29. , 29.25, 29.5 , 29.75, 30. , 30.25, 30.5 , 30.75, 31. , 31.25,
31.5 , 31.75, 32. , 32.25, 32.5 , 32.75, 33. , 33.25, 33.5 , 33.75,
34. , 34.25, 34.5 , 34.75, 35. , 35.25, 35.5 , 35.75, 36. , 36.25,
36.5 , 36.75, 37. , 37.25, 37.5 , 37.75, 38. , 38.25, 38.5 ])array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
array([[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]],
[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]],
[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]],
[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]])gpm2009 = awash_da1
gpm2019 = awash_da2
imd2009 = ds3_ind_sm.RAINFALL
imd2019 = ds4_ind_sm.RAINFALL
# using interp_like
imd2009_interp = imd2009.interp_like(gpm2009)
imd2019_interp = imd2019.interp_like(gpm2019)
Error
# 2009 and 2019 GPM variation comparison to IMD
err2009 = gpm2009 - imd2009_interp
err2019 = gpm2019 - imd2019_interp
err2009.sel(season = 'MAM').plot()
<matplotlib.collections.QuadMesh at 0x7fd780967c10>
err2019.sel(season = 'MAM').plot()
<matplotlib.collections.QuadMesh at 0x7fd7807bca00>
# Plotting overall error for 2009 and 2019
# Plotting 2009 and 2019 RMSE
fig = plt.figure(figsize=(20, 15))
fig.tight_layout()
ax = fig.add_subplot(1, 2, 1, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
err2009.mean(dim='season').plot(cbar_kwargs = {"orientation":"horizontal", "pad":0.05})
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("2009", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
ax = fig.add_subplot(1, 2, 2, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
err2019.mean(dim='season').plot(cbar_kwargs = {"orientation":"horizontal", "pad":0.05})
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("2019", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
fig.suptitle('GPM data error compared to IMD gridded data (in mm)', fontsize=20, y=0.78)
plt.savefig('./images/err.png')
/home/aditya/.local/share/virtualenvs/atms_python-xEvIgfwt/lib/python3.9/site-packages/dask/array/numpy_compat.py:39: RuntimeWarning: invalid value encountered in true_divide x = np.divide(x1, x2, out) /home/aditya/.local/share/virtualenvs/atms_python-xEvIgfwt/lib/python3.9/site-packages/dask/array/numpy_compat.py:39: RuntimeWarning: invalid value encountered in true_divide x = np.divide(x1, x2, out)
# Plotting 2009 seasonal error
fig = plt.figure(figsize=(20, 15))
fig.tight_layout()
titles = ["MAM", "JJA", "SON", "DJF"]
for i,season in enumerate(titles):
ax = fig.add_subplot(2, 2, i+1, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
err2009.sel(season=titles[i]).plot()
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("Months"+ " " + "-" + " " + "("+titles[i]+")", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
fig.suptitle('GPM data error compared to IMD gridded data (in mm)-2009', fontsize=20, y=0.95)
plt.savefig('./images/err2009.png')
# Plotting 2019 seasonal error
fig = plt.figure(figsize=(20, 15))
fig.tight_layout()
titles = ["MAM", "JJA", "SON", "DJF"]
for i,season in enumerate(titles):
ax = fig.add_subplot(2, 2, i+1, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
err2019.sel(season=titles[i]).plot()
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("Months"+ " " + "-" + " " + "("+titles[i]+")", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
fig.suptitle('GPM data error compared to IMD gridded data (in mm)-2019', fontsize=20, y=0.95)
plt.savefig('./images/err2019.png')
RMSE
# Resetting Colormap selection for RMSE
xr.set_options(cmap_divergent='bwr', cmap_sequential='CMRmap') # divergent doesn't matter here
<xarray.core.options.set_options at 0x7f846d800ac0>
rmse2009 = np.sqrt((err2009 * err2009).mean(dim = 'season'))
rmse2019 = np.sqrt((err2019 * err2019).mean(dim = 'season'))
# Plotting 2009 and 2019 RMSE
fig = plt.figure(figsize=(20, 15))
fig.tight_layout()
ax = fig.add_subplot(1, 2, 1, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
rmse2009.plot(cbar_kwargs = {"orientation":"horizontal", "pad":0.05})
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("2009", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
ax = fig.add_subplot(1, 2, 2, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
rmse2019.plot(cbar_kwargs = {"orientation":"horizontal", "pad":0.05})
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("2019", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
fig.suptitle('GPM RMSE compared to IMD gridded data (in mm)', fontsize=20, y=0.78)
plt.savefig('./images/rmse.png')
# define a function to calculate rmse error for given dataarray error value
def rmse_calc(da_err, season):
"""
The RMSE calc function calculates the rmse from given input error value
and also takes a season string as input for selecting the seasonal
mean whose rmse needs to be calculated
"""
months = ['MAM','JJA','SON','DJF']
if season == 'MAM':
rmse = np.sqrt((da_err * da_err).sel(season = months[0]))
elif season == 'JJA':
rmse = np.sqrt((da_err * da_err).sel(season = months[1]))
elif season == 'SON':
rmse = np.sqrt((da_err * da_err).sel(season = months[2]))
elif season == 'DJF':
rmse = np.sqrt((da_err * da_err).sel(season = months[3]))
else:
print("ERROR : Please enter a correct season value")
return rmse
# Plotting rmse error for seasonal means for 2009
fig = plt.figure(figsize=(20, 15))
fig.tight_layout()
titles = ["MAM", "JJA", "SON", "DJF"]
for i,season in enumerate(titles):
ax = fig.add_subplot(2, 2, i+1, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
rmse_calc(err2009, titles[i]).plot()
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("Months"+ " " + "-" + " " + "("+titles[i]+")", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
fig.suptitle('2009 GPM RMSE compared to IMD gridded data (in mm)', fontsize=20, y=0.95)
plt.savefig('./images/rmse2009.png')
# Plotting rmse error for seasonal means for 2019
fig = plt.figure(figsize=(20, 15))
fig.tight_layout()
titles = ["MAM", "JJA", "SON", "DJF"]
for i,season in enumerate(titles):
ax = fig.add_subplot(2, 2, i+1, projection=ccrs.PlateCarree())
ax.set_extent([67, 98, 7, 36], crs=ccrs.PlateCarree())
rmse_calc(err2019, titles[i]).plot()
gridliner = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
gridliner.top_labels = False
gridliner.bottom_labels = True
gridliner.left_labels = True
gridliner.right_labels = False
gridliner.ylines = False # you need False
gridliner.xlines = False # you need False
ax.set_title("Months"+ " " + "-" + " " + "("+titles[i]+")", pad=10, fontsize=14)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
fig.suptitle('2019 GPM RMSE compared to IMD gridded data (in mm)', fontsize=20, y=0.95)
plt.savefig('./images/rmse2019.png')
# Function to convert the IMD daily data to one month data
# we do the same with GPM for consistency
def daily_to_monthly(da):
temp = da.groupby('time.month').mean(dim='time')
times = pd.date_range(start = "2009-01-15", end = "2010-01-15",freq='M')
output = xr.DataArray(
temp,
coords={
"time": times,
"lon": temp.lon,
"lat": temp.lat
},
dims=["time", "lat", "lon"],
)
return output
# get the time series dataarrays using the above function. Here "ts" is timeseries
gpm2009_ts = daily_to_monthly(ds1_ind.precipitation)
imd2009_ts = daily_to_monthly(ds3_ind.RAINFALL)
gpm2019_ts = daily_to_monthly(ds2_ind.precipitation)
imd2019_ts = daily_to_monthly(ds4_ind.RAINFALL)
# Time series 2009
fig= plt.figure(figsize=(9,6))
imd2009_ts.mean(dim='lat').mean(dim='lon').plot.line("-*",color='orange', label = 'IMD')
gpm2009_ts.mean(dim='lat').mean(dim='lon').plot.line("b-^", label = 'GPM')
plt.legend()
plt.title('2009 precipitaion time-series for India', fontdict={"size":18})
plt.ylabel('Precipitation (in mm)', fontdict={"size":14})
plt.xlabel('time',fontdict={"size":14})
plt.savefig('./images/time_series2009.png')
# Time series 2019
fig = plt.figure(figsize=(9,6))
imd2019_ts.mean(dim='lat').mean(dim='lon').plot.line("-*",color='orange', label = 'IMD')
gpm2019_ts.mean(dim='lat').mean(dim='lon').plot.line("b-^", label = 'GPM')
plt.legend()
plt.title('2019 precipitaion time-series for India', fontdict={"size":18})
plt.ylabel('Precipitation (in mm)', fontdict={"size":14})
plt.xlabel('time',fontdict={"size":14})
plt.savefig('./images/time_series2019.png')